***************
* Agriculture *
***************

/* Legend:
	1 = thousand hectares
	2 = thousand tonnes
	3 = tonnes per hectare
*/

* Import agricultural data
	foreach a in rice maize soybean wheat {
	foreach x in 1 2 3 {
		clear all
		import excel "orig/OECD/agriculture.xlsx", sheet("`a'`x'") firstrow
		destring cou, replace
		save "data/`a'`x'.dta", replace
	}
	}

* Merge with main data set
	clear all
	use "data/final_data_set.dta"
	xtset id year
	foreach a in rice maize soybean wheat {
	foreach x in 1 2 3 {
		merge m:m cou year using "data/`a'`x'.dta", keep(match master)
		drop _merge
		erase "data/`a'`x'.dta"
	}
	}
	xtset id year
	
* Generate scaled output measures
	foreach a in rice maize soybean wheat {
	foreach x in 1 2 3 {
		gen `a'`x'_rgdpna = `a'`x' / rgdpna
	}
	}

* Generate aggregate output measures
	foreach x in 1 2 3 {
		gen agr`x' = rice`x' + maize`x' + soybean`x' + wheat`x'
		gen agr`x'_rgdpna = ln(agr`x' / rgdpna)
		gen agr`x'_1000 = agr`x' / 1000 
	}
	
* Summary statistics	
	// Panel A. Full sample
		tabstat agr1_rgdpna agr2_rgdpna agr3_rgdpna agr1_1000 agr2_1000 agr3_1000, ///
		columns(statistics) statistics(mean sd p25 p50 p75) format(%9.4f)
	// Panel B. G7 countries
		tabstat agr1_rgdpna agr2_rgdpna agr3_rgdpna agr1_1000 agr2_1000 agr3_1000 if dum_G7 == 1, ///
		columns(statistics) statistics(mean sd p25 p50 p75) format(%9.4f)		
	// Panel C. Non-G7 countries
		tabstat agr1_rgdpna agr2_rgdpna agr3_rgdpna agr1_1000 agr2_1000 agr3_1000 if dum_G7 == 0, ///
		columns(statistics) statistics(mean sd p25 p50 p75) format(%9.4f)		
	
* FE regressions
	// relative measures
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_rgdpna dum_G7, robust fe(year) cluster(id)
				eststo
		}
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_rgdpna dum_G7, robust fe(year) cluster(year)
				eststo
		}
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear
	// absolute measures
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_1000 dum_G7, robust fe(year) cluster(id)
				eststo
		}
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_1000 dum_G7, robust fe(year) cluster(year)
				eststo
		}
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear

* Merge with weather data
	merge m:m cou year using "data/weather.dta"
	drop _merge

	duplicates drop id year, force
	xtset id year

	// weather variables: DP01 DP10 DP1X PRCP DP01_i DP10_i DP1X_i PRCP_i

	// generate precipitation variable that excludes extreme rain days
		gen PRCP1 = PRCP - DP1X

	// generate rates of growth and dummies
		foreach x in DP01 DP10 DP1X PRCP PRCP1 DP01_i DP10_i DP1X_i PRCP_i {
			gen g_`x' = ln(`x')-ln(l.`x')
			gen dum_`x' = (d.`x' < 0)
			* replace dum_`x' = . if d.`x' == . // dubious command, creates issues below (missing obs)
			replace dum_`x' = . if pd == .
		}

	// generate alternative dummy using moving average of past three years	
		foreach x in PRCP PRCP1 {
		rangestat (mean) `x', interval(year -4 -1) by(id)
		gen dum_`x'_mean = (`x' - `x'_mean < 0)
		replace dum_`x' = . if pd == .
		}
		// correlations across dummies
			pwcorr dum_P*, sig	
		
	// exogeneity assumption
		pwcorr rtfpna f_rtfpna zcci dum_PRCP dum_PRCP1 DP1X PRCP PRCP1 dum_DP1X, sig
		xtreg rtfpna zcci dum_PRCP1 i.year if pd != ., robust fe // choice: dum_PRCP dum_PRCP1 DP1X PRCP PRCP1 dum_DP1X 
		reg rtfpna zcci dum_PRCP1
		
	// auxiliary regression
		xtreg cci dum_PRCP i.year if pd != ., robust fe // choice: dum_PRCP dum_PRCP1 DP1X PRCP PRCP1 dum_DP1X 
		predict cci_w, xb
		predict cci_nw, r
		* NOTE 1: the precipitation dummy allows for comparisons across countries as it abstracts from levels
		* NOTE 2: to use three-year moving average, replace with "dum_PRCP_mean"

	// standardized
		foreach x in cci_w cci_nw {
		egen `x'mean = mean(`x')
		egen `x'sd = sd(`x')
		gen z`x' = (`x'-`x'mean) / `x'sd
		drop `x'sd `x'mean
		}

	xtset id year
	generate newid = id
	tsset newid year
		
* Regressions with weather data
	// relative measures
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_rgdpna dum_PRCP, robust fe(year) cluster(id)
				eststo
		}
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_rgdpna dum_PRCP, robust fe(year) cluster(year)
				eststo
		}
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear
	// absolute measures
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_1000 dum_PRCP, robust fe(year) cluster(id)
				eststo
		}
		foreach x in 1 2 3 {
			qui: regxfe agr`x'_1000 dum_PRCP, robust fe(year) cluster(year)
				eststo
		}
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear
		
* Interactions between weather and G7		
		gen dum_PRCP_dum_G7 = dum_PRCP * dum_G7
	// relative measures
		qui: regxfe agr3_rgdpna dum_PRCP dum_G7 dum_PRCP_dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear
	// absolute measures
		qui: regxfe agr3_1000 dum_PRCP dum_G7 dum_PRCP_dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear		
		
* Correlations
	pwcorr agr3 agr3_rgdpna dum_PRCP dum_G7, sig
	
* Final table
	qui: regxfe agr1_rgdpna dum_PRCP, robust fe(year) cluster(id)
		eststo
	qui: regxfe agr1_rgdpna dum_G7, robust fe(year) cluster(id)
		eststo
	qui: regxfe agr1_rgdpna dum_PRCP dum_G7, robust fe(year) cluster(id)
		eststo
	qui: regxfe agr1_rgdpna dum_PRCP dum_G7 dum_PRCP_dum_G7, robust fe(year) cluster(id)
		eststo
	esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
	t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
	drop()
	estimates clear
	
	tabulate country if dum_PRCP != . & agr1_rgdpna != .	
		********* Australia, Canada, New Zealand, Switzerland, United Kingdom, United States
	tabstat year if dum_PRCP != . & agr1_rgdpna != ., statistics(min max)
		********* sample starts in 1990
			

	
***********************
* Missed days at work *
***********************

* Import data on missed days of work
	clear all
	import excel "orig/OECD/missed_days_at_work.xlsx", sheet("`a'`x'") firstrow
	destring cou, replace
	save "data/missed_days_at_work.dta", replace

* Merge with main data set
	clear all
	use "data/final_data_set.dta"
	xtset id year
	merge m:m cou year using "data/missed_days_at_work.dta", keep(match master)
	drop _merge
	erase "data/missed_days_at_work.dta"
	duplicates drop id year, force
	xtset id year
	
* Generate scaled output measures
	gen missed_rgdpna = ln(missed / rgdpna)

* Summary statistics	
	// Panel A. Full sample
		tabstat missed missed_rgdpna, ///
		columns(statistics) statistics(mean sd p25 p50 p75) format(%9.4f)
	// Panel B. G7 countries
		tabstat missed missed_rgdpna if dum_G7 == 1, ///
		columns(statistics) statistics(mean sd p25 p50 p75) format(%9.4f)		
	// Panel C. Non-G7 countries
		tabstat missed missed_rgdpna if dum_G7 == 0, ///
		columns(statistics) statistics(mean sd p25 p50 p75) format(%9.4f)		
	
* FE regressions
	// relative measures
		qui: regxfe missed_rgdpna dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear
	// absolute measures
		qui: regxfe missed dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear

* Merge with weather data
	merge m:m cou year using "data/weather.dta"
	drop _merge

	duplicates drop id year, force
	xtset id year

	// weather variables: DP01 DP10 DP1X PRCP DP01_i DP10_i DP1X_i PRCP_i

	// generate precipitation variable that excludes extreme rain days
		gen PRCP1 = PRCP - DP1X

	// generate rates of growth and dummies
		foreach x in DP01 DP10 DP1X PRCP PRCP1 DP01_i DP10_i DP1X_i PRCP_i {
			gen g_`x' = ln(`x')-ln(l.`x')
			gen dum_`x' = (d.`x' < 0)
			* replace dum_`x' = . if d.`x' == . // dubious command, creates issues below (missing obs)
			replace dum_`x' = . if pd == .
		}

	// generate alternative dummy using moving average of past three years	
		foreach x in PRCP PRCP1 {
		rangestat (mean) `x', interval(year -4 -1) by(id)
		gen dum_`x'_mean = (`x' - `x'_mean < 0)
		replace dum_`x' = . if pd == .
		}
		// correlations across dummies
			pwcorr dum_P*, sig	
		
	// exogeneity assumption
		pwcorr rtfpna f_rtfpna zcci dum_PRCP dum_PRCP1 DP1X PRCP PRCP1 dum_DP1X, sig
		xtreg rtfpna zcci dum_PRCP1 i.year if pd != ., robust fe // choice: dum_PRCP dum_PRCP1 DP1X PRCP PRCP1 dum_DP1X 
		reg rtfpna zcci dum_PRCP1
		
	// auxiliary regression
		xtreg cci dum_PRCP i.year if pd != ., robust fe // choice: dum_PRCP dum_PRCP1 DP1X PRCP PRCP1 dum_DP1X 
		predict cci_w, xb
		predict cci_nw, r
		* NOTE 1: the precipitation dummy allows for comparisons across countries as it abstracts from levels
		* NOTE 2: to use three-year moving average, replace with "dum_PRCP_mean"

	// standardized
		foreach x in cci_w cci_nw {
		egen `x'mean = mean(`x')
		egen `x'sd = sd(`x')
		gen z`x' = (`x'-`x'mean) / `x'sd
		drop `x'sd `x'mean
		}

	xtset id year
	generate newid = id
	tsset newid year
		
* Regressions with weather data
	// relative measures
		qui: regxfe missed_rgdpna dum_PRCP dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear
	// absolute measures
		qui: regxfe missed dum_PRCP dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear

* Interactions between weather and G7		
		gen dum_PRCP_dum_G7 = dum_PRCP * dum_G7
	// relative measures
		qui: regxfe missed_rgdpna dum_PRCP dum_G7 dum_PRCP_dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear
	// absolute measures
		qui: regxfe missed dum_PRCP dum_G7 dum_PRCP_dum_G7, robust fe(year) cluster(id)
			eststo
		esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
		t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
		drop()
		estimates clear		
		
* Simple correlations		
	pwcorr missed dum_PRCP dum_G7, sig
	
* Final table
	qui: regxfe missed_rgdpna dum_PRCP, robust fe(year) cluster(id)
		eststo
	qui: regxfe missed_rgdpna dum_G7, robust fe(year) cluster(id)
		eststo
	qui: regxfe missed_rgdpna dum_PRCP dum_G7, robust fe(year) cluster(id)
		eststo
	qui: regxfe missed_rgdpna dum_PRCP dum_G7 dum_PRCP_dum_G7, robust fe(year) cluster(id)
		eststo
	esttab, cells(b(star fmt(%9.4f)) t(fmt(%9.2f))) style(fixed) title(Table_15_c) ///
	t starlevels(* .10 ** .05 *** .01) b(4) scalars(r2o)     ///
	drop()
	estimates clear
	
	tabulate country if dum_PRCP != . & missed_rgdpna != .
		********* all except New Zealand
	tabstat year if dum_PRCP != . & missed_rgdpna != ., statistics(min max)
		********* sample starts in 2010
		